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A SIMPLE LINEAR RESPONSE CLOSURE APPROXIMATION FOR SLOW 
DYNAMICS OF A MULTISCALE SYSTEM WITH LINEAR COUPLING 

RAFAIL V. ABRAMOV 

Abstract. Many applications of contemporary science involve multiscale dy- 
namics, which are typically characterized by the time and space scale separation 
of patterns of motion, with fewer slowly evolving variables and much larger set 
^ I of faster evolving variables. This time-space scale separation causes direct numer- 

ical simulation of the evolution of the dynamics to be computationally expensive, 

00 I due both to the large number of variables and the necessity to choose a small 

discretization time step in order to resolve the fast components of dynamics. In 
this work we propose a simple method of determining the closed model for slow 

j^f^ • variables alone, which requires only a single computation of appropriate statistics 

^^ I for the fast dynamics with a certain fixed state of the slow variables. The method 

is based on the first-order Taylor expansion of the averaged coupling term with 
respect to the slow variables, which can be computed using the linear fluctuation- 
C^ . dissipation theorem. We show that, with simple linear coupling in both slow and 

fast variables, this method produces quite comparable statistics to what is exhib- 
ited by a complete two-scale model. The main advantage of the method is that it 
applies even when the statistics of the full multiscale model cannot be simulated 
^ i due to computational complexity, which makes it practical for real-world large 

00 ' scale applications. 
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1. Introduction 



Multiscale dynamics are common in applications of contemporary science, such 
as geophysical science and climate change prediction dUlHUlllHl]]. Multiscale 
dynamics are typically characterized by the time and space scale separation of 
patterns of motion, with (typically) fewer slowly evolving variables and much 
5^ \ larger set of faster evolving variables. This time-space scale separation causes 

direct numerical simulation of the evolution of the dynamics be computationally 
expensive, due both to the large number of variables and the necessity to choose a 
small discretization time step in order to resolve the fast components of dynamics. 

In the climate change science the situation is further complicated by the fact 
that climate is characterized by the long-term statistics of the slow variables, which, 
under small changes of parameters (such as the solar radiation forcing, green- 
house gas content, etc) change over even longer time scale than the motion of 
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the slow variables themselves. In this situation, where long-term statistics of the 
slow motion patterns need to be captured, the direct forward time integration of 
the most comprehensive global circulation models (GCM) is subject to enormous 
computational expense. 

As a more computationally feasible alternative to direct forward time inte- 
gration of the complete multiscale model, it has long been recognized that, if 
a closed simplified model for the slow variables alone is available, one could 
use this closed slow-variable model instead to simulate the statistics of the slow 
variables. Numerous closure schemes were developed for multiscale dynamical 
systems iri0l[T3ll20l - |23l , which are all based on the averaging principle over the 
fast variables Il25ll32ll33ll . Some of the methods (such as those in H20l - |23l ) replace 
the fast nonlinear dynamics with suitable stochastic processes IMI or conditional 
Markov chains [10], while others [13] provide direct closure by suitable tabula- 
tion and curve fitting. However, it seems that all these approaches require either 
extensive computations to produce a closed model (for example, [10,13] require 
multiple simulations of fast variables alone with different fixed states of slow 
variables), or somewhat ad hoc determination of closure coefficients by matching 
areas under the time correlation functions [20-23]. 

In this work we propose a simple method of determining the closed model for 
slow variables alone, which requires only a single computation of appropriate 
statistics for the fast dynamics with a certain fixed state of the slow variables. 
The method is based on the first-order Taylor expansion of the averaged coupling 
term with respect to the slow variables, which, as we show, can be computed 
using the linear fluctuation-dissipation theorem in]-l8l [T9ll26| . We show through 
the computations with the appropriately rescaled two-scale Lorenz 96 model [HI 
that, with simple linear coupling in both slow and fast variables, this method 
produces quite comparable statistics to what is exhibited by the slow variables 
of the complete two-scale Lorenz model. The main advantage of the method is 
its simplicity and easiness of implementation, partly due to the fact that the fast 
dynamics need not be explicitly known (that is, the fast dynamics can be provided 
as a "black-box" algorithm), and the parameters of the closed model for the slow 
variables are determined from the appropriate statistics of the fast variables for a 
given fixed state of the slow variables. Additionally, the method can be applied 
even when the statistics for slow variables of the full multiscale model are not 
available due to computational expense. 

The manuscript is structured as follows. In Section |2] we outline the theoretical 
grounds for the algorithm. In Section |3] we describe the two-scale Lorenz model 
Ill3[ll7[[18l , appropriately rescaled so that the means and variances of both the fast 
and slow variables are near zero and one, respectively [4]. Section H] contains the 
results of numerical experiments with both the two-scale Lorenz model and the 
reduced set of equations for slow variables only, comparing various statistics of 
the time series. Section |5] summarizes the results of this work. 
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2. The linear response closure approximation 
Consider a two-scale system of differential equations of the form 

(2.1) ^ = ^(^'^)' 1J = ^(^'^)' 

where x = x{t) G R^^ are the slow variables, y = y{t) G R^y are the fast vari- 
ables, and F and G are Nx and Ny vector-valued functions of x and y, respectively. 
Here and below, we assume that the fast variables y are "unresolved", that is, the 
computation of the dynamics for y requires such a small time discretization step, 
that the direct computation of (|2.1[) for long time intervals is practically infeasible. 
We also assume that the dynamics for y are fast enough for the system in (|2.1|) to 
be approximated by the averaged dynamics for x, given by 

(2.2) ^ = (F)(x), (F)(x) = |^^^F(x,z)d^,(z), 

for finite times (for a more detailed description of the averaging formalism, see 
Pll4l l25ll32ll33| ). Above, }ix is the invariant probability measure of the limiting fast 
dynamics given by 

dz 

(2.3) — = G(x,z), 

dT 

where the solution is given by the flow z(t) = ^Jzq, while x is a fixed constant 
parameter for (|2.3|) and, consequently, }ix- Below we assume that (F)(x) varies 
smoothly with respect to x (it is known that for stochastically driven systems this 
property is generic, and for deterministic dynamics this happens when }ix is an 
SRB measure ||12ll27l - |29l l35]). Under the ergodicity assumption for }ix, one can 
replace the measure average with time average: 

(2.4) (F)(x) = lim - f Fix, zi T))dT, 

r^co r Jo 

where z(t) is a long-term trajectory of (|2.3[) . 

While suitable multiscale methods exist (see | [TT| and references therein) where 
the computation of the average in (|2.4[) is performed relatively rarely, even that 



might not be computationally feasible in complex multiscale systems with many 
fast variables. For the purpose of this work, here we assume that the computation 
of (|2.3|) is practically feasible only for a single choice of the constant parameter 
X = X*, where x* is a suitable point, in the vicinity of which the motion occurs, 
such as the mean state of the original dynamics in (|2.1|) , or a nearby state. A poor 
man's approach in this case is to compute the approximate average with respect 
to Hx*, which is a zero order approximation: 

(2.5) {F){x)= I F{x,z)d}ix*{z) + 0{\\x-x*\\). 

Jwry 
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Here, one has to compute the time average in (|2.4|) only once, for the time series 
of (|2.3|) corresponding to x = x*. However, as recently found in HI], this approx- 



imation may fail to capture the chaotic properties of the slow variables in (|2.1|) 
Here we propose the following first order correction: 



(2.6) 






+ 0(\\x- 



,*l|2^ 



where }i'^* is the derivative of fix with respect to x, computed at the point x = x*. 
It is shown in BH] that the average of F with respect to }i'^* can be computed as the 
corresponding infinite time linear response to infinitesimal perturbation of x (for 
details on the linear response theory, see l|TI48l [T9lE6ll and references therein): 

^^•^^ L^^ ^^""'^^ "^^^^^^^ = L L^y a^(^''/'x*2)n*,, g^(x*,z) d^,*(z) ds, 

where dF/dy is the Jacobian with respect to the second argument of F, and Tj* ^ 
is the tangent map of (|2.3|) : 



(2.8) 






dz 



Under the assumption of ergodicity of fix*' (|2.7|| can be written as the integral 
of the time autocorrelation function, which is computed along a single long-time 



trajectory z*(t) of (|23l) with x = x*: 



(2.9) 



/ F(x,z)duL 



dG 



L ,'S"J„ 3r(^.3-(T+»))n-M37(-*.^-w)dT 



ds. 



which is also done only once along the same long-term trajectory of (|2.3|) , as for 
the time average in (|2.5|) (for practical purposes, it can be assumed that the time 
autocorrelation function decays sufficiently fast to replace the improper integral 
from zero to infinity with a proper integral up to a sufficiently long time). While 
formally valid, the above formula can be unsuitable for practical computation 
due to the fact that usually the limiting fast dynamics in (|2.3[) are strongly chaotic 
and mixing with large positive Lyapunov exponents. The presence of large pos- 
itive Lyapunov exponents causes numerical instability in the computation of the 
tangent map T of (|2.3|) for long response times, and, as a result, the infinite-time 
linear response in (|2.9|) cannot be computed. Below we consider the special setting 
for (|2.1[) with linear coupling (which is common in geophysical sciences), and use 
the quasi-Gaussian linear response approximation [[l-:3,.5-8,.19J for the practical 
computation of (|2.9[) . 
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2.1. Special case with linear coupling. Here we consider the special setting of 
(|2.1|) with linear coupling between x and y: 

(2.10) F{x, y) = fix) + Lyy, G{x, y) = g{y) + L,x, 

where / and g are nonlinear vector functions of x and y, respectively, and Lx and 
Ly are constant matrices of suitable sizes. For this simplified setting, observe that 
dG/dx = Lx, dF/dy = Ly, and, therefore, the approximate averaged system is 
given by 

c\x 

(2.11a) — = fix) + Lyz* + LyR*Lx{x - x*), 

(2.11b) z*= lim- r z*(r)dr, 

r^oo r Jo 

(2.11c) R^ 



lim - / ^^w^-^dT ds, 

where z{t) is the solution of the fast limiting system 

dz 

(2.12) ^ = g{z) + Lxx 

with X specified as a constant parameter, and z*{t) corresponds to x = x*. The 
formula in (|2.11c|) is usually unsuitable for direct numerical computation due to 
rapidly growing Lyapunov exponents at fast scales. Instead, one can use the 
quasi-Gaussian linear response approximation, where (|2.11c|) becomes the inte- 
gral of the time autocorrelation function under the assumption of the Gaussian 
invariant measure for (|2.12|) . With linear coupling, the measure-averaged linear 
response formula (|2.7[) becomes 

(2.13) ^^^ F{x, z) d^;* (z) = Ly (j^ ^^^ Tl.^, dfi^^ (z) ds j Lx- 
Now, the Gaussian probability density is given by 

(2.14) pliz) = (271)-^ detL*-i exp ("^(^ " z-*)^E*-i(z - z*)) , 
where E* is the long-time covariance matrix of (|2.12|) , computed at the point x*: 

(2.15) Z* = lim - r(z*(T) - z*)(z*(t) - z*)^dT. 

r->oo r Jo 

Recalling that T^* ^ = dcp^^z/dz and replacing d}ix*{z) = p^(z)dz, we obtain, 
after integration by parts over z, 

(2.16) j^^^F{x,z)dM^) = -h (/f /^,^</>^*z^dzds) Lx. 
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Computing the derivative of (|2.14[) with respect to z, we obtain 



(2.17) 1^^^ I{x,z) Ay!A^) = Ly {^j^ j^^^ f,.z{z - z*fph{z) dzds) 1.*-'L,. 



Switching back to time-averaging, we obtain 



(2.18) 



R% 



F{x,z) d}i'x*{z] 



.y r [lim- rz*(T + s)(z*(T)-Z*)^dT 
^ Jo r^co r Jo 



/O i^r^-00 r Jo 

and, therefore, R* in (|2.11a|) can now be computed as 



dsV 



(2.19) R* = lim- / z*(t + s)(z 

Jo r^oo r Jo 



Z*)^dT 



dsL' 



For details, see (IHIHSHSHTII). 

For even better precision of the linear response computation, one can also use 
the blended linear response approximation |I5]-^, however, in this work we do 
not implement it, as it is shown that for the model and regimes considered, the 
quasi-Gaussian approximation is already quite precise. 

Even with the linear coupling, the function z{x) (the dependence of the mean 



state of (|2.12[) on x) is not generally linear. Thus, the validity of the linear approx- 
imation in (|2.11a|) depends on the influence (or lack thereof) of the nonlinearity of 
the function z{x). While rigorous estimates of the validity of the linear approxi- 
mation in (|2.11a|) can hardly be provided in general case, here, instead, we try to 
justify it by comparing the fast limiting system in (|2.12|) to the Ornstein-Uhlenbeck 
process ||3T1 . Consider an Ornstein-Uhlenbeck process of the form 



(2.20) 



dz 
dr 



-r- = — r(z — m) + LxX + (T 



dWr 

"dT' 



where m is a constant Ny-vector, F is a constant Ny x Ny positive-definite matrix, 
Wf is a i<C-dimensional Wiener process, cr is a constant Ny x K matrix, and x is, 
as in (|2.12|| , is a constant parameter. Then, it is easy to see that the difference 
between the statistical mean states of (|2.20|) corresponding to x and x* is 

(2.21) zou-Zou = ^~^^x{x-x*), 

which is valid for {x — x*) of an arbitrary norm. At the same time, by the re- 
gression theorem II26L the time correlation function of p.20|) with x = x* is given 

by 



(2.22) 



{zhu{t + s){zhuit) 






f) = exp(-sr)L 



OU' 



where Eqjj is the covariance matrix of the Ornstein-Uhlenbeck process in (|2.20|) 



for X = X*. Thus, according to (|2.19)) , the infinite-time linear response operator 
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for (|2.20|) is computed as 



(2.23) Rqu = / exp(-sr) ds 



By comparing (|2.23[) with (|2.21[) , one can see that, for the Ornstein-Uhlenbeck pro- 



cess, the quasi-Gaussian linear response formula in (|2.19|) is exact for an arbitrarily 
large perturbation {x — x*). Hence, if the nonlinear process in (|2.12|) behaves sta- 
tistically similarly to the Ornstein-Uhlenbeck process in (|2.20[) , the averaged sys- 
tem in (|2.11a|l can be expected to behave statistically similarly to the slow part of 
(|2.10|) . Below we numerically test the approximation for slow variables in the spe- 
cial setting with linear coupling using the two-scale Lorenz model I.1..2l141 [T01[T31i17i1 . 



3. The two-scale Lorenz model 

Here we choose the two-scale forced damped Lorenz model 
for the computational study of the dynamical properties of a two-scale slow-fast 
process with generic features of climate-weather systems, such as the presence 
of linearly unstable waves, strong nonlinearity, forcing, dissipation, chaos and 
mixing. The two-scale forced damped Lorenz model is given by 

(3.1a) Xi = Xi_i{xi+i - Xi-2) -Xi + Fx- ^Yl yi,j, 

1 A 

(3.1b) 1/g- = - [y/,y+i(yfj_i - yfj+2) - yi,j + Fy] + -^x,-, 

where 1 <i < Nx, 1 < / < /• The following notations are adopted above: 

• X is the set of the slow variables of size Nx- The following periodic bound- 
ary conditions hold for x: Xj+j^x = ^u 

• y is the set of the fast variables of size Ny = Nx} where / is a positive 
integer. The following boundary conditions hold for y: yi+N^,i = yfj arid 

yi,i+} = yi+i,j' 

• Fx and Fy are the constant forcing parameters; 

• \x arid Ay are the coupling parameters; 

• e is the time scale separation parameter. 

Originally in ||10||13[|T7| there was no constant forcing F term in the equation for x- 
variables in (|3.1|) , however, in its absence the behavior of the y-variables is strongly 
dissipative [1,2J. Here, as in [2J, we add a constant forcing Fy in the right-hand 
side of the second equation in (|3.1|) to induce the strongly chaotic behavior of the 
y-variables with large positive Lyapunov exponents. 
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3.1. Rescaled Lorenz model. In the Lorenz model (|3.1|) , Fx and Fy regulate the 
chaos and mixing of the x and y variables, respectively ||ll|2l,|4l. However, the 
mean state and mean energy are also affected by the changes in forcing, which 
affects the mean and energy trends in coupling for the fixed coupling parameters. 
To adjust the effect of coupling independently of forcing, here we rescale the 
Lorenz model as in Ifl9| . Consider the uncoupled Lorenz model 

(3.2) —Xi = x,_i(x,+i - Xi_2) -Xi + F 

with the same periodic boundary conditions as above flE\. Observe that the long 



term statistical mean state x and the standard deviation /3 in (|3.2[) are the same for 



all Xi due to the translational invariance. Now, we rescale x and t as 

T 

(3.3) Xi = x + ^Xi, t = -, 

where the new variables x have zero mean state and unit standard deviation. In 
the rescaled variables, the Lorenz model becomes 

d 1 _ F — X 

(3.4) —Xi = X;-l(Xf+l - Xi-2) + -o [x{Xi+l - Xi-2) - Xi] + —P2-' 

where x and /3 are, of course, the functions of F. In addition to setting the mean 
state and variance of X; to zero and one, respectively, due to the time rescaling 
the autocorrelation functions of z acquire roughly identical time scaling for any F 
(for details, see O). Here, we similarly rescale the two-scale Lorenz model from 

1 F — X A ^ 
(3.5a) Xi = X/_i(x/+i - X/_2) + — (x(x,+i - Xi-2) - Xi) + -^ ^ J] i//,y. 



yr,i 



(3.5b) 



yi,j+iiyi,i-i - y 1,1+2) + ^ (y(y/,i-i - y 1,1+2) - yi,j) + 

Py 



PI 



~r X;, 

e 



where x, y, /3x and /Sy are the long term means and standard deviations of the cor- 
responding uncoupled system in (|3.2|) with either F^ or Fy set as a constant forcing. 
In the rescaled Lorenz model (|3.5|) , the values of Fx and Fy do not significantly af- 
fect the mean state, mean energy and the time scale for both the slow variables x 
and fast variables y, and mostly regulate the mixing and non-Gaussianity of the 
probability distributions of the long-term time series. 
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The Lorenz model (|3.5|) represents the setting with linear coupling as in (|2.10|) , 
with f,g,Lx and Ly given by 

1 F — X 

(3.6a) fi{x) = Xi-i{xi+i - Xi-2) + — (x(xf+i - Xi-2) - Xi) + "" ^ ' 



8t,i{y) = - 
(3.6b) 



1 

y/,i+i(y/j-i - y/,7+2) + ^ (y(yz,;-i - y/,7+2) - yzj) + 

Py 



/3^ 



y 



(3.6c) ^^ = T^^' ^y = -^^' hm = ^tj- 

Now, one can immediately see that for the Lorenz model in (|3.5)) , the limiting 
system in (|2.3)) is given by 



(3.7) ^-^(^) + T^'- 

while the approximate averaged system around the point x* is given by 

(3.8) ^ = fix) - ^Lz* - ^LR^L^ix - x*). 

Above, the parameter e in the denominator suggests that the first-order correction 
term is of order e~^, which is misleading, because R* scales proportionally to 
e (since it is the integral of the time autocorrelation function with rate of decay 
proportional to £~^). For simplicity of computation, for the Lorenz model we 
rescale the fast limiting dynamics by e to bring the time scale to order 1 since the 
scaling parameter e is known explicitlyli which yields 

dz 

(3.9) — = eg(z)+A,L^x, 

for the fast limiting dynamics for fixed x, while for the averaged dynamics we 
obtain 

(3.10) ^ = fix) - ^Lz* - ^LR*L^ix - x*), 

where R* is computed from the time series of p.9|) using the quasi-Gaussian 
approximation in (|2.19|) . 



Generally, when no scaling parameter is available explicitly in l|2.1l l, one can still multiply the 
right-hand side of I l2.3t by a heuristic small parameter to bring the time scale of Il2.3t to order 
1, since the averaged dynamics are invariant with respect to the £-scaling of the limiting fast 
dynamics (for details, see [4J). Or, equivalently, use appropriately small time step and averaging 
window for 
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4. Numerical experiments 

Here we present a numerical study of the proposed approximation for slow dy- 
namics, applied to the rescaled Lorenz model in (|3.5|| . We compare the statistical 
properties of the slow variables for the three following systems: 

(1) The complete rescaled Lorenz system from (|3.5|) ; 

(2) The approximation for slow dynamics alone from (|3.10|) : 



(3) The poor man's version of (|3.10|) with the first-order correction term R* set 
to zero (further referred to as the "zero-order" system). 
The fixed parameter x* for the computation of R* was set to the long-term mean 
state X of (|3.5|) (in practical situations, a rough estimate could be used). The quasi- 
Gaussian approximation in (|2.19|) is used to compute the first order correction 
term R*. While it is practically impossible to compute the improper integrals 
in (|2.19|) for infinite upper limit, in practice we use sufficiently long (but finite) 



limits of integration. In particular, for all computational results presented below, 
the correction term R* is computed numerically as 

(4.1) R* = ^ r™"' r"\*{T + s){z*{T)-z*fdTdsZ*-\ 

i-av Jo Jo 

where the averaging time window Tav equals 10000 time units, while the cor- 
relation time window Tcorr equals 50 time units (it was observed that the time 
autocorrelation function in (|4.1|) decays essentially to zero within the 50 time-unit 
window for all studied regimes). The mean state z* and the covariance matrix E* 
are also computed by time-averaging with the same averaging window of 10000 
time units. 

Due to translational invariance of the studied models, the statistics are invariant 
with respect to the index shift for the variables x,. For diagnostics, we monitor 
the following long-term statistical quantities of Xf. 

a. The probability density functions (PDF), computed by bin-counting. A PDF 
gives the most complete information about the one-point statistics of x„ as it 
shows the statistical distribution of Xj in the phase space. 

b. The time autocorrelation functions {xi{t)xi{t + s)), where the time average is 
over t, normalized by the variance (xf) (so that it always starts with 1). 

c. The time cross-correlation functions (xf(f)xf_|_i(f + s)), also normalized by the 



variance (xj) 



d. The energy autocorrelation function 

{xj(t)xj{t + s)) 



m 



This energy autocorrelation function measures the non-Gaussianity of the pro- 
cess (it is identically 1 for all s if the process is Gaussian, such as the Ornstein- 
Uhlenbeck process). For details, see [2] 
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The success (or failure) of the proposed approximation of the slow dynamics 
depends on several factors. First, as the quasi-Gaussian linear response formula 
is used for the computation of R*, the precision will be affected by the non- 
Gaussianity of the fast dynamics. Second, it depends how linearly the mean state 
z for the fast variables depends on the slow variables x. Here we observe the 
limitations of the proposed approximation by studying a variety of dynamical 
regimes of the rescaled Lorenz model in (|3.5|) . The following dynamical regimes 
are studied: 

• Nx = 20, } = A (so that Ny = 80). Thus, the number of the fast variables is 
four times greater than the number of the slow variables. 

• e = 0.01. The time scale separation of two orders of magnitude is consistent 
with typical real-world geophysical processes (for example, the annual and 
diurnal cycles of the Earth's atmosphere). 

• Ax = Ay = 0.3,0.4. These values of coupling are chosen so that they 
are neither too weak, nor too strong (although 0.3 is weaker, and 0.4 is 
stronger). Recall that the standard deviations of both x, and i/; , variables 
are approximately 1, and, thus, the contribution to the right-hand side 
from coupled variables is weaker than the self-contribution, but still of the 
same order. 

• Fx = 6, 16. The slow forcing Fx adjusts the chaos and mixing properties 
of the slow variables, and in this work it is set to a weakly chaotic regime 
Fx = 6, and strongly chaotic regime Fx = 16. 

• Fy = 8, 12. The fast forcing adjusts the chaos and mixing properties of the 
fast variables. Here the value of Fy is chosen so that the fast variables are 
either moderately chaotic for Fy = 8, or more strongly chaotic for Fy = 12. 

In Figure |4H we show the probability density functions and time autocorrelation 
functions for the limiting fast dynamics in (|3.9|) , with the parameters Nx = 20, 
Ny = 80, Fx = 6, Ax = Ay = 0.3, and two values of the fast forcing: f y = 8 and 
Fy = 12. Observe that the PDFs are not Gaussian (although close to it), and have 
nonzero skewness. The time autocorrelation functions decay slower for Fy = 8 
and faster for Fy = 12, indicating slower and faster mixing, respectively. In other 
regimes, these PDFs and autocorrelation functions look very similar to what is 
presented in Figure 14. 1[ 

Another point we would like to emphasize before presenting the results of the 
computational study, is that the matrices R* and, consequently, LRLJ in (|3.10|) 
are not diagonal. In Figure |42] we display both R* and LRLJ for the dynamical 
regime with Fx = 6, Fy = 8, and A^ = Ay = 0.4 (only the central columns 
of R* and LRLJ are displayed with diagonal elements corresponding to zero 
horizontal coordinates of the plots, as both matrices are translation-invariant). 
Observe that there are significant off-diagonal entries in both matrices. Both R* 
and LRLJ are positive-definite in the presented regime (the lowest eigenvalues of 
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Nx=20, N =80, Fx=6, ;ix=0.3, ?Ly=0.3 



Nx=20, Ny=80, Fx=6, ;i>,=0.3, ;iy=0.3 

Fy=12 - 




Figure 4.1. The probability density functions and time autocorrela- 
tion functions of the limiting fast dynamics in (|3.9[) for Fx = 6, and 
Ax = Ay = 0.3. 
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Figure 4.2. The matrices R* (left) and LR*L^ (right) of the reduced 
dynamics in (|3.10|) for Fx = 6, Fy = 8, and Ax = Ay = 0.4. Only 
a single column of each matrix is displayed with its diagonal ele- 
ment corresponding to zero horizontal coordinate, as the matrices 
are translation-invariant. 



their symmetric parts are 6.314 • 10 ~^ and 0.8814, respectively), and, as a result, the 
linear correction term in (|3.10|) causes damping effect on the reduced dynamics 
(for more details about chaotic properties of reduced dynamics, see [4J). For other 
regimes, the matrices are similar to the presented regime (that is, substantial off- 
diagonal entries are present), and we do not display those here. 
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Figure 4.3. PDFs, A^ = A^ = 0.3. 



4.1. Probability density functions of the slow dynamics. In Figures l43l and l4!4l 
we show the probability density functions of the slow dynamics for the full two- 
scale Lorenz model, the reduced closed model for the slow variables alone in 
(|3.10|) , and its poor man's zero order version without the linear correction term. 
Observe that for the more weakly coupled regimes with A^ = Ay = 0.3 the PDFs 
look rather similar, however, it can be seen that the reduced model with the correc- 
tion term reproduces the PDFs much closer to those of the full two-scale Lorenz 
model, than the zero-order model. In the more strongly coupled regime with 
Ax = Ay = 0.4 the situation tilts even more in favor of the reduced model with 
linear correction term in (|3.10|) : observe that for the weakly chaotic regime with 
Fx = 6 the PDFs of the full two-scale Lorenz model have three sharp peaks, in- 
dicating strong non-Gaussianity The reduced model in (|3.10|) reproduces these 
peaks, while its zero-order version fails. In addition, in Table 1411 we show the 
L2-errors in PDFs between the full two-scale Lorenz model and the two reduced 
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Figure 4.4. PDFs, A^ = Ay = 0.4. 

models. Observe that, generally, the reduced system with the linear correction 
term in (|3.10|) produces more precise results than its poor man's version without 
the correction term. 



4.2. Time autocorrelation functions of the slow dynamics. In Figures 14.51 and 
14.61 we show the time autocorrelation functions of the slow dynamics for the full 
two-scale Lorenz model, the reduced closed model for the slow variables alone in 
(|3.10|) , and its poor man's zero order version without the linear correction term. 
Observe that for the more weakly coupled regimes with A^ = Ay = 0.3 the time 
autocorrelation functions look similar, yet the reduced model with the correction 
term reproduces the time autocorrelation functions more precisely than the zero- 
order model. In the more strongly coupled regime with A^ = Ay = 0.4 the 
difference between the reduced model in (|3.10|) and its poor man's zero-order 
version is even more drastic: observe that for the weakly chaotic regime with 
Fx = 6 the time autocorrelation functions of the full two-scale Lorenz model do 
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Ax,v = 0-3/ Fy = 8 


\x,y = 0.3, Fy = 12 




Red. 
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= 16 
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1.469 - 10-2 


Fx 
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Red. 
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Red. 


Z.O. 


F, 
Fx- 
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= 16 


0.1022 
3.725 • 10-3 


8.857-10-2 
2.703 - 10-2 


Fx 
Fx- 


= 6 
= 16 


9.28 - 10-2 
5.885 - 10-3 


0.1113 
3.209 - 10-2 



Table 4.1. L2-errors between the PDFs of the slow variables of the 
full two-scale Lorenz model and the two reduced models. Nota- 
tions: "Red." stands for "Reduced" (that is, (I3l0ll ), and "Z.O." 
stands for "Zero-order", the poor man's version of (|3.10|) . 



\x,y = 0.3, Fy = 8 


\x,y = 0.3, Fy = 12 




Av 


Z.O. 




Av 


Z.O. 


Fx 
Fx-- 


= 6 
= 16 


5.841 - 10-2 
4.079 - 10-2 


0.1211 
5.342 - 10-2 


Fx- 


= 6 
= 16 


6.539 - 10-2 
1.559 - 10-2 


0.1572 
7.396 - 10-2 




\x,y = 0.4, f^ = 8 


\x,y = 0.4, Fy = : 


L2 




Av 


Z.O. 




Av 


Z.O. 


Fx = 6 
Fx = l6 


5.538 - 10-2 
8.534 - 10-2 


0.3677 
0.1355 


Fx 
Fx- 


= 6 
= 16 


0.2981 
4.835 - 10-2 


0.3986 
0.1482 



Table 4.2. L2-errors between the time autocorrelation functions of 
the slow variables of the full two-scale Lorenz model and the two 
reduced models. Notations: "Red." stands for "Reduced" (that is, 
(|3.10|) ), and "Z.O." stands for "Zero-order", the poor man's version 
of (IsiUl) . 



not exhibit decay (indicating very weak mixing), and the reduced model in (|3.10|) 
reproduces the autocorrelation functions of the full two-scale Lorenz model rather 
well, while its zero-order version fails. In addition, in Table |4^ we show the L2- 
errors in time autocorrelation functions (for the correlation time interval of 20 
time units, as in Figures 14.51 and 14. 6|) between the full two-scale Lorenz model and 
the two reduced models. Observe that, generally, the reduced system with the 
linear correction term in (|3.10|) produces more precise results than its poor man's 
version without the correction term. 



4.3. Time cross-correlation functions. In Figures 14.71 and 14.81 we show the time 
cross-correlation functions of the slow dynamics for the full two-scale Lorenz 
model, the reduced closed model for the slow variables alone in (|3.10)) , and its 
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Figure 4.5. Time autocorrelation functions, A^ = Ay = 0.3. 



poor man's zero order version without the linear correction term. Observe that for 
the more weakly coupled regimes with A^ = Ay = 0.3 the time cross-correlation 
functions look similar, however, it is seen that the reduced model with the cor- 
rection term reproduces the time cross-correlation functions more precisely than 
the zero-order model. In the more strongly coupled regime with A^ = Ay = 0.4, 



the reduced model in (|3.10|) becomes much more precise than its poor man's 
zero-order version: here, for the weakly chaotic regime with F^ = 6 the time 
cross-correlation functions of the full two-scale Lorenz model do not exhibit de- 
cay (indicating very weak mixing), and the reduced model in (|3.10|) reproduces the 
cross-correlation functions of the full two-scale Lorenz model rather well, while 
its zero-order version fails. In addition, in Table 14.31 we show the L2-errors in 
time cross-correlation functions (for the correlation time interval of 20 time units, 
as in Figures 14.71 and 14.8)1 between the full two-scale Lorenz model and the two 
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Figure 4.6. Time autocorrelation functions, A^ = Ay = 0.4. 



reduced models. Observe that, generally, the reduced system with the linear cor- 
rection term in (|3.10|) produces more precise results than its poor man's version 
without the correction term. 



4.4. Energy autocorrelation functions. In Figures 14.91 and 14.101 we show the en- 
ergy autocorrelation functions of the slow dynamics for the full two-scale Lorenz 
model, the reduced closed model for the slow variables alone in (|3.10|) , and its 
poor man's zero order version without the linear correction term. Observe that 
for the more weakly coupled regimes with A^ = Ay = 0.3 the energy autocor- 
relation functions look similar, although the reduced model with the correction 
term reproduces the energy autocorrelation functions more precisely than the 
zero-order model. In the more strongly coupled regime with A^ = Ay = 0.4, 



the reduced model in (|3.10|) is much more precise than its poor man's zero-order 
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Figure 4.7. Time cross-correlation functions, Xx = )^v = 0.3. 



version: here, for the weakly chaotic regime with F^ = 6 the energy autocorre- 
lation functions of the full two-scale Lorenz model is significantly sub-Gaussian, 
and the reduced model in (|3.10|) reproduces the sub-Gaussianity of the energy 
autocorrelation functions of the full two-scale Lorenz model rather well, while its 
zero-order version fails. In addition, in Table 14.41 we show the L2-errors in en- 
ergy autocorrelation functions (for the correlation time interval of 20 time units, 
as in Figures 14.91 and |4.10|) between the full two-scale Lorenz model and the two 
reduced models. Observe that, generally, the reduced system with the linear cor- 
rection term in (|3.10|) produces more precise results than its poor man's version 
without the correction term. 
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Figure 4.8. Time cross-correlation functions. Ax = '^v = 0.4. 



5. Summary 

In this work we develop a simple method of constructing the closed reduced 
model for slow variables of a multiscale model with linear coupling, which re- 
quires only a single computation of the mean state and the time autocorrelation 
function for the fast dynamics with a fixed state of the slow variables, which is 
located in the region where the slow dynamics evolve (here, the mean state of the 
slow dynamics is used). The method is based on the first-order Taylor expansion 
of the averaged coupling term with respect to the slow variables, which is com- 
puted using the linear fluctuation-dissipation theorem. We demonstrate through 
the computations with the appropriately rescaled two-scale Lorenz 96 model [4| 
that, with simple linear coupling in both slow and fast variables, the developed 
reduced model produces quite comparable statistics to what is exhibited by the 
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Table 4.3. L2-errors between the time cross-correlation functions of 
the slow variables of the full two-scale Lorenz model and the two 
reduced models. Notations: "Red." stands for "Reduced" (that is, 
(|3.10|) ), and "Z.O." stands for "Zero-order", the poor man's version 
of (IsiUl) . 
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Table 4.4. L2-errors between the energy autocorrelation functions 
of the slow variables of the full two-scale Lorenz model and the two 
reduced models. Notations: "Red." stands for "Reduced" (that is, 
(|3.10|) ), and "Z.O." stands for "Zero-order", the poor man's version 
of (l3igil . 



complete two-scale Lorenz model. Below we outline the main advantages of the 
new method: 

• The reduced model is simple. It requires only the mean state, covariance 
and the time autocorrelation function for the fast variables, computed for a 
single fixed state of the slow variables. Since only the statistics of the time 
series of the fast dynamics are needed, the structure of the right-hand side 
of the equations for the fast variables need not be known explicitly - this 
part of dynamics can be provided as a "black box". Also, the structure of 
the nonlinear x-dependent part of the right-hand side of the equations for 
the slow variables need not be known to construct the approximation, and 
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Figure 4.9. Energy autocorrelation functions, A^ = Ay = 0.3. 



existing computational routines can be used for it. The only correction in 
the forward time-stepping routine is the linear correction term. 

• The reduced model is a priori. It lacks parameters which have to be ad- 
justed a posteriori to "fit" the statistical properties of the full multiscale 
dynamics. In fact, statistical properties of the full multiscale dynamics 
need not be known to construct the reduced model (although certain sta- 
tistics of the fast variables with an appropriate fixed slow state need to be 
computed). 

• The reduced model is parsimonious. It requires only a simple linear correc- 
tion to achieve consistently better performance than that of a correspond- 
ing zero-order model. 
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Figure 4.10. Energy autocorrelation functions, A^ = Ay = 0.4. 



• The reduced model is practical. It can be implemented even when the 
statistics for the slow variables of the complete multiscale dynamics can- 
not be obtained due to its computational complexity (although a rough 
estimate of the mean state, or a nearby state is needed), which makes the 
approach potentially suitable for comprehensive global circulation models 
in geophysics. Additionally, existing zero-order models (such as the T21 
barotropic model LZnl4„30J) can be retrofitted with the linear correction 
term. 

In the future work, the author intends to collaborate with geophysicists to create 
more realistic reduced models for geophysical dynamics, including retrofitting 
existing closed models for slow dynamics with the linear correction term. 
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